home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Source Code / Libraries / DCLAP 4j / SeqPups / apps / clustalw.src / readmat.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-12-17  |  7.3 KB  |  331 lines  |  [TEXT/R*ch]
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include <stdlib.h>
  4. #include <string.h>
  5. #include <ctype.h>
  6. #include "clustalw.h"
  7. #include "matrices.h"
  8.  
  9. /*
  10.  *   Prototypes
  11.  */
  12.  
  13. void     init_matrix(void);
  14. int     get_matrix(int *matptr, int *xref,
  15.                      int matrix[NUMRES][NUMRES], int neg_flag);
  16. int     read_user_matrix(char *filename, int *usermat, int *xref);
  17. int     getargs(char *inline1,char *args[],int max);
  18.  
  19. /*
  20.  *   Global variables
  21.  */
  22.  
  23. extern char     *amino_acid_codes;
  24. extern int     gap_pos1, gap_pos2;
  25. extern int     max_aa;
  26. extern int     def_aa_xref[], aa_xref[], pw_aa_xref[];
  27. extern int     usermat[], pw_usermat[];
  28. extern int     mat_avscore;
  29. extern int     debug;
  30.  
  31. void init_matrix(void)
  32. {
  33.  
  34.    char c1,c2;
  35.    int i, j, maxres;
  36.  
  37.    max_aa = strlen(amino_acid_codes)-2;
  38.    gap_pos1 = max_aa + 1;          /* code for gaps inserted by clustalw */
  39.    gap_pos2 = gap_pos1 + 1;           /* code for gaps already in alignment */
  40.  
  41. /*
  42.    set up cross-reference for default matrices hard-coded in matrices.h
  43. */
  44.    for (i=0;i<NUMRES;i++) def_aa_xref[i] = -1;
  45.  
  46.    maxres = 0;
  47.    for (i=0;c1=amino_acid_order[i];i++)
  48.      {
  49.          for (j=0;c2=amino_acid_codes[j];j++)
  50.           {
  51.            if (c1 == c2)
  52.                {
  53.                   def_aa_xref[i] = j;
  54.                   maxres++;
  55.                   break;
  56.                }
  57.           }
  58.          if ((def_aa_xref[i] == -1) && (amino_acid_order[i] != '*'))
  59.             {
  60.                 fprintf(stderr,
  61.                 "Warning: residue %c in matrices.h is not recognised\n",
  62.                                        amino_acid_order[i]);
  63.             }
  64.      }
  65.  
  66. }
  67.  
  68. int get_matrix(int *matptr, int *xref, int matrix[NUMRES][NUMRES], int neg_flag)
  69. {
  70.    double f;
  71.    int gg_score = 1;
  72.    int gr_score = 0;
  73.    int i, j, k, ix = 0;
  74.    int ti, tj;
  75.    int  maxres, min;
  76.    int av1,av2,av3, max;
  77.  
  78. /*
  79.    default - set all scores to 0
  80. */
  81.    for (i=0;i<=max_aa;i++)
  82.       for (j=0;j<=max_aa;j++)
  83.           matrix[i][j] = 0;
  84.  
  85.    ix = 0;
  86.    maxres = 0;
  87.    for (i=0;i<=max_aa;i++)
  88.     {
  89.       ti = xref[i];
  90.       for (j=0;j<=i;j++)
  91.        {
  92.           tj = xref[j]; 
  93.           if ((ti != -1) && (tj != -1))
  94.             {
  95.                k = matptr[ix];
  96.                if (ti==tj)
  97.                   {
  98.                      matrix[ti][ti] = k * 100;
  99.                      maxres++;
  100.                   }
  101.                else
  102.                   {
  103.                      matrix[ti][tj] = k * 100;
  104.                      matrix[tj][ti] = k * 100;
  105.                   }
  106.                ix++;
  107.             }
  108.        }
  109.     }
  110.  
  111.    av1 = av2 = av3 = 0;
  112.    for (i=0;i<=max_aa;i++)
  113.     {
  114.       for (j=0;j<=i;j++)
  115.        {
  116.            av1 += matrix[i][j];
  117.            if (i==j)
  118.               {
  119.                  av2 += matrix[i][j];
  120.               }
  121.            else
  122.               {
  123.                  av3 += matrix[i][j];
  124.               }
  125.        }
  126.     }
  127.  
  128.    --maxres;
  129.    av1 /= (maxres*maxres)/2;
  130.    av2 /= maxres;
  131.    av3 /= ((float)(maxres*maxres-maxres))/2;
  132. if (debug>1) fprintf(stderr,"average mismatch score %d\n",(pint)av3);
  133. if (debug>1) fprintf(stderr,"average match score %d\n",(pint)av2);
  134. if (debug>1) fprintf(stderr,"average score %d\n",(pint)av1);
  135.  
  136.   min = max = matrix[0][0];
  137.   for (i=0;i<=max_aa;i++)
  138.     for (j=1;j<=i;j++)
  139.       {
  140.         if (matrix[i][j] < min) min = matrix[i][j];
  141.         if (matrix[i][j] > max) max = matrix[i][j];
  142.       }
  143. /*
  144.    if requested, make a positive matrix - add -(lowest score) to every entry
  145. */
  146.   if (neg_flag == FALSE)
  147.    {
  148.  
  149. if (debug>1) fprintf(stderr,"min %d max %d\n",(pint)min,(pint)max);
  150.       if (min < 0)
  151.         {
  152.            for (i=0;i<=max_aa;i++)
  153.             {
  154.               ti = xref[i];
  155.               if (ti != -1)
  156.                 {
  157.                  for (j=0;j<=max_aa;j++)
  158.                    {
  159.                     tj = xref[j];
  160. /*
  161.                     if (tj != -1) matrix[ti][tj] -= (2*av3);
  162. */
  163.                     if (tj != -1) matrix[ti][tj] -= min;
  164.                    }
  165.                 }
  166.             }
  167.         }
  168. /*
  169.        gr_score = av3;
  170.        gg_score = -av3;
  171. */
  172.  
  173.    }
  174.  
  175.   mat_avscore = -av3;
  176.  
  177.  
  178.   for (i=0;i<gap_pos1;i++)
  179.    {
  180.       matrix[i][gap_pos1] = gr_score;
  181.       matrix[gap_pos1][i] = gr_score;
  182.       matrix[i][gap_pos2] = gr_score;
  183.       matrix[gap_pos2][i] = gr_score;
  184.    }
  185.   matrix[gap_pos1][gap_pos1] = gg_score;
  186.   matrix[gap_pos2][gap_pos2] = gg_score;
  187.   matrix[gap_pos2][gap_pos1] = gg_score;
  188.   matrix[gap_pos2][gap_pos1] = gg_score;
  189.  
  190.   maxres += 2;
  191.  
  192.   return(maxres);
  193. }
  194.  
  195.  
  196. int read_user_matrix(char *filename, int *usermat, int *xref)
  197. {
  198.    double f;
  199.    int gg_score = 0;
  200.    int  numargs;
  201.    int i, j, k, ix1, ix = 0;
  202.    FILE *fd;
  203.    char codes[NUMRES];
  204.    char inline1[1024];
  205.    char *args[NUMRES+4];
  206.    char c1,c2;
  207.    int  maxres, min;
  208.    int sum;
  209.  
  210.    if (filename[0] == '\0')
  211.      {
  212.         fprintf(stderr,"Error: comparison matrix not specified\n");
  213.      }
  214.    else
  215.      {
  216.         if ((fd=fopen(filename,"r"))==NULL) 
  217.           {
  218.              fprintf(stderr,"Error: cannot open %s\n", filename);
  219.              return(0);
  220.           }
  221.  
  222.         maxres = 0;
  223.         while (fgets(inline1,1024,fd) != NULL)
  224.           {
  225.              if ((inline1[0] == '\0') || (inline1[0] == '#')) continue;
  226. /*
  227.    read residue characters.
  228. */
  229.              k = 0;
  230.              for (j=0;j<strlen(inline1);j++)
  231.                {
  232.                   if (isalpha((int)inline1[j])) codes[k++] = inline1[j];
  233.                   if (k>NUMRES)
  234.                      {
  235.                         fprintf(stderr,"Error: too many entries in %s\n",filename);
  236.                         return(0);
  237.                      }
  238.                }
  239.              codes[k] = '\0';
  240.              break;
  241.           }
  242.  
  243.         if (k == 0) 
  244.           {
  245.              fprintf(stderr,"Error: wrong format in %s\n",filename);
  246.              return(0);
  247.           }
  248.  
  249. /*
  250.    cross-reference the residues
  251. */
  252.    for (i=0;i<NUMRES;i++) xref[i] = -1;
  253.  
  254.    maxres = 0;
  255.    for (i=0;c1=codes[i];i++)
  256.      {
  257.          for (j=0;c2=amino_acid_codes[j];j++)
  258.            if (c1 == c2)
  259.                {
  260.                   xref[i] = j;
  261.                   maxres++;
  262.                   break;
  263.                }
  264.          if ((xref[i] == -1) && (codes[i] != '*'))
  265.             {
  266.                 fprintf(stderr,"Warning: residue %c in %s not recognised\n",
  267.                                        codes[i],filename);
  268.             }
  269.      }
  270.  
  271.  
  272. /*
  273.    get the weights
  274. */
  275.  
  276.         ix = ix1 = 0;
  277.         while (fgets(inline1,1024,fd) != NULL)
  278.           {
  279.              if (inline1[0] == '\n') continue;
  280.              numargs = getargs(inline1, args, maxres+1);
  281.              if (numargs == 0)
  282.                {
  283.                   fprintf(stderr,"Error: wrong format in %s\n", filename);
  284.                   return(0);
  285.                }
  286.              for (i=0;i<=ix;i++)
  287.                {
  288.                   if (xref[i] != -1)
  289.                     {
  290.                        f = atof(args[i]);
  291.                        usermat[ix1++] = (int)(f);
  292.                     }
  293.                }
  294.              ix++;
  295.           }
  296.         if (ix != maxres+1)
  297.           {
  298.              fprintf(stderr,"Error: wrong format in %s\n", filename);
  299.              return(0);
  300.           }
  301.      }
  302.  
  303.     fclose(fd);
  304.  
  305.   maxres += 2;
  306.  
  307.   return(maxres);
  308. }
  309.  
  310. int getargs(char *inline1,char *args[],int max)
  311. {
  312.  
  313.     char    *inptr;
  314. #ifndef MAC
  315.     char    *strtok(char *s1, const char *s2);
  316. #endif
  317.     int    i;
  318.  
  319.     inptr=inline1;
  320.     for (i=0;i<=max;i++)
  321.     {
  322.         if ((args[i]=strtok(inptr," \t\n"))==NULL)
  323.             break;
  324.         inptr=NULL;
  325.     }
  326.     if (i!=max) return(0);
  327.  
  328.     return(i);
  329. }
  330.  
  331.